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A CONTRACTING- INTERVAL PROGRAM FOR THE 


DANILEWSKI METHOD* 

By James D. Harris 
Langley Research Center 

SUMMARY 

The concept of contracting- interval programs is applied to the Danilewski method 
for finding the eigenvalues of a matrix. The development is a three-step process in 
which (1) a contracting- interval program is developed for the reduction of a matrix to 
Hessenberg form, (2) a contracting- interval program is developed for the reduction of a 
Hessenberg matrix to colleague form, and (3) the characteristic polynomial with interval 
coefficients is readily obtained from the interval of colleague matrices, and this interval 
polynomial is then factored into quadratic factors so that the eigenvalues may be 
obtained. 

To develop a contracting-interval program for factoring this polynomial with inter- 
val coefficients it is necessary to have an iteration method which converges even in the 
presence of controlled rounding errors. 

A theorem is stated giving sufficient conditions for the convergence of Newton's 
method when both the function and its Jacobian cannot be evaluated exactly but when the 
errors can be made proportional to the square of the norm of the difference between the 
previous two iterates. This theorem is applied to prove the convergence of the gener- 
alization of the Newton- Bairstow method that is used to obtain quadratic factors of the 
characteristic polynomial. 


This research was originally presented to the Faculty of the School of Engi- 
neering and Applied Science, University of Virginia, in partial fulfillment of the 
requirements for the degree of Doctor of Philosophy (Applied Mathematics). 
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SYMBOLS 


Standard matrix notation is used throughout; that is, elements are designated by 
lower case letters corresponding to the matrix symbol and with appropriate subscripts. 

A real n by n matrix 

A A after row and column interchanges 

B matrix containing rounding errors made in reduction of H to F 

E matrix containing rounding errors made in converting A to H 

F colleague matrix in chapter in 

F defined as NM* - I in chapter II 

f(x) vector whose components are functions of x 

f|(x) approximation to f 

H upper Hessenberg matrix 

I identity matrix 

Ij^, identity matrix with columns j and j' interchanged 

J Jacobian of f 

J j. approximation to J 

M inverse of N 

N lower triangular matrix 

n order of matrix A 

P(x) polynomial of degree n 

Pj(x) polynomials satisfying recurrence relation 
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Pj real numbers satisfying recurrence relation 

R large positive number 

n-dimensional real- coordinate space 
u(s,t),v(s,t) functions defined following algorithm H of chapter IV 
V upper triangular matrix 

W inverse of V 

superdiagonal element of colleague matrix 
diagonal element in colleague matrix 
6f(x) positive continuous function of x 

6J(x) positive continuous function of x 

6P(x) polynomial of degree n - 1 with positive coefficients 

X constant with same dimensions as 1/x 

Notation used is as follows: 


where A is a matrix of order n, ||a|| = m^ ^ ^ 




A ^ B if A and B are vectors or matrices, this means that each component 

of A is less than or equal to corresponding component of B 

[a ± 6A] {^A^Ia - 6A = A^ ^ A + 6A} or interval with midpoint A and half- 
width 6A 


f(S) 


asterisk following a symbol denotes approximation; for instance, u*(s,t) 
is an approximation to u(s,t) 

{f(X)|x e S} where f is a function and S c R*^ 
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n 

. 2 . 

[P ± 6P] 


sup 


l|x|| 

{x|s} 


greatest integer less than 


n 

2 


where P(x) = + a^^ jPj^ + . . . + &qP^ and 

6P = Sa^^ + 5aj^_j^Pj + . . . + (Sa^ g 0 for each i), 

denotes set of all polynomials p'(x) = a' + . . . + alp , + a, 

n 1 n-1 I 

that a| € [a^ ± Sa^J (i = l,2,...,n) 
product symbol 


supremum (least upper bound) 


where, x is a vector 


Xi 


X2 


X = max ( X, 


ISi^n V ^ 


set of all X such that statement s is true 


min-^|s^ least element of indicated set 
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,P such 
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CHAPTER I 


INTRODUCTION 

This paper presents the results of an investigation of techniques for applyir^ the 
concept of contracting-interval programs (ref. 1) to the Danilewski method for finding the 
eigenvalues of a matrix. 


Justification 

Most data used for input to a computer program are measured data or are the 
result of truncating exact data. These inexact data are then processed by a computer 
program which makes rounding errors. Thus, much analysis may be required to deter- 
mine how closely the computer solution approximates the exact solution of the problem. 

Contracting- interval programs present an alternative to this procedure. The user 
of a contracting-interval program can supply an interval which is known to contain the 
exact data. The contracting- interval program then gives him an interval, each element 
of which is the exact answer corresponding to some element in the input interval. This 
is as good a result as he could hope to obtain. 

Notice that contracting- interval programs are not the same as expanding- interval 
programs. An expanding- interval program computes an output interval which contains 
all of the answers corresponding to elements of the input interval. However, there may 
be elements in the output interval of an expanding- interval program which do not corre- 
spond to any element in the input interval. For a further discussion of the differences 
between contracting- and expanding- interval programs, see reference 1. 

The Danilewski method is particularly appropriate for testing the concept of 
contracting- interval programs, since it is a method which is very sensitive to rounding 
errors. When executed by using a fixed-precision arithmetic, it frequently gives unreli- 
able results because of the rounding errors, but a contracting- interval program based on 
the Danilewski method yields exact results. Also, the Danilewski method is a fast 
method for obtaining eigenvalues and, even with the increased precision that a 
contracting- interval program requires, it may prove to be competitive in speed with 
other more inherently stable methods (ref. 2). 

The Danilewski method used in the present paper is not the classical Danilewski 
method but is a three- step process whereby (1) a matrix is reduced to Hessenberg form 
by a similarity transformation, (2) the resulting Hessenberg matrix is reduced to col- 
league form (ref. 3) by a similarity transformation, and (3) the characteristic polynomial 
is obtained from the colleague matrix and the roots of this polynomial are calculated. 



As Wilkinson (ref. 4) has shown, it is in this second step that problems develop. The 
resulting colleague matrix may be significantly more ill-conditioned than was the corre- 
sponding Hessenberg matrix. 

Previously, contractir^- interval programs have been developed for solving linear 
equations (ref. 1) and for finding real roots of a polynomial (ref. 5), but a contracting- 
interval program has not been attempted for a problem so complex as the matrix- 
eigenvalue problem. 


Preview of Remaining Chapters 

Any exact computation defines a mapping f from one vector space to 

another vector space R™. If S c r”, then f(S) = {f(X)|X e S} . If [A ± 5A] C 
then f[A ± 6A] is not necessarily an interval. However, a contracting-interval pro- 
gram computes an interval [B ± 6B] such that [B ± 6B] c f [A ± 6A]. 


A two-step process is used here to develop contracting-interval programs. The 
first step obtains B such that B = f(A') where I A - A'| = and R is a large 
positive number. During this step the program must have the ability to control the 
rounding errors made in each arithmetic operation. The second step finds 6B such 


that [B ± 6B] c f 
|A' - A"| S ^ 1 


A' ± ^ ~ ^ 6A 


R 


; that is, if B’ e [B ± 6B], then B’ 


R 


6A. 


f(A") where 


Contracting- interval programs have the property (ref. 1) that, given contracting- 
interval programs for f and g, there is automatically obtained a contracting- interval 
program for the composition of f and g. 


In chapter II a contracting- interval program is developed for reducing an 
interval [A ± 6A] of real matrices to an interval [H ± 6H] of upper Hessenberg 
matrices by a similarity transformation. Most of the material in chapter II is new, 
although the algorithm for obtaining it is based on the one given in reference 4 on 
page 357, 

In chapter IH a con tractir^- interval program for reducing an interval [H ± 6H] 
of upper Hessenberg matrices to an interval [F ± 6F] of colleague matrices by a simi- 
larity transformation is developed. The concept of a colleague matrix that is used here 
is a generalization of the one given in reference 3. The form of the colleague matrix as 
shown in chapter III was suggested by B. A. Chartres. Much of the work in chapter IH 
is based on work done in reference 6 althoi^h, in reference 6, Chartres was concerned 
with obtaining a companion matrix F which was similar to some H’ e [H ± 6H]. 

From [F ± 6F] a set of polynomials [P ± 6P] is obtained, such that each element is 
the characteristic polynomial of some A' e [A ± 6A]. 
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In chapter IV an algorithm is developed for obtaining an interval [bj ± 6b J and 
the intervals jsj^ ± 6s^ and jt^ ± 6t^ ^i = 1 , 2 ,..., ^ j, where is the greatest 

integer in ^ and n is the degree of P such that, if s' e [sj^ ± 6s^, t' e [tj ± 6t^, 
and b’ e [bj ± 6b then 


^ao + bjPj) jy ^Pg - SjPj - tj^ e [P ± 6P] 


If n is even, then bj = 6bj^ = 0. The method used is a generalization of the Newton- 
Bairstow method. In appendix A, a proof is given that the algorithm developed in 
chapter IV leads to convergence. 

From chapter IV it may be seen that Pg - sjPj - tj is a factor of the character- 
istic polynomial of some matrix A' e [A ± 6A] for sj e ± 6siJ and t^ e [t^ ± 6t^ 

^i = 1,2,...,^!^ j. Thus, the roots of this quadratic polynomial are eigenvalues of the 
matrix A*. Intervals about the eigenvalues are not determined in this study. 
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CHAPTER n 


REDUCTION OF AN INTERVAL OF REAL MATRICES TO AN INTERVAL 
OF UPPER HESSENBERG MATRICES 

In chapter II an algorithm for converting an interval of real matrices into an 
interval of upper Hessenberg matrices by a similarity transformation is described. 

Let A be a real matrix of interval midpoints, and let 6A be a matrix whose 
elements are the interval half-widths. Let H denote an upper Hessenberg matrix; 
that is, H has the form 



Let N denote a matrix of the form 



Let 6H be an upper Hessenberg matrix of all positive elements. The matrices H 
and 6H are determined such that if H' e [H ± 6H], then H' is similar to some 
A' e [A ± 6A]. To do this, an N of the form given above is found such that if 
H' e [H ± 6H], then H’ = N'^A’N for some A' e [A ± 6A]; that is, 

[H ± 6H] C N“1[A ± 6A]N. 
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Since A is the midpoint of the original interval, it is desirable that the midpoint 
of the interval of Hessenberg matrices satisfy H = N"^AN for some N; that is, 

NH = AN (1) 

Such an H can be calculated by using exact arithmetic in the following algorithm as 
long as no zero elements are generated along the subdiagonal of H (ref. 4, p. 358): 


Algorithm A 
For r = 1,2, ...,n 



Now some provision must be made so that the algorithm applies in the case when 
an r ~ ® calculated. This is done in the following manner: 

A lgorithm B 
For r = 1,2, ...,r 

For i = 1,2, ...,r 

n i-1 

*^ir = ^ir + ^ ^ik'^kr ' 2/ "ik\r 

k=r+l k=2 
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If r < n, then 


For i = r+l,r+2,...,n 


CT. = a- + 
1 ir 


k=r+l 


Find j such that 


a. = max ^ cr. 


nr 1 r+2|’-’ I 

Interchange 
If 0, then 

For i = r +l,r+2,...,n 
Interchange 

For i = 1,2, ■■.,n 
Interchange 

For i = l,^,.,.,r 
Interchange 

^r+l,r “ °^r+l 

For i = r+2,r+3,...,n 


i,r+l h 


r+l,r 


If ^r+1 

For i = r+2,r+3,...,n 


\r+l=0 


With algorithm B, an N and an H are obtained, but equation (1) is no longer 
satisfied. 


Effect of Interchanges 


The effect of interchanges is shown in theorem 1 in which 1^., denotes the iden- 
tity matrix with columns i and j' interchanged. When Ijj, multiplies a matrix on 
the left, the result is to interchange rows j and j' in that matrix. When multi- 

plies a matrix on the right, the result is to interchange columns j and j' in that 
matrix. 

Theorem 1 — Given an n x n matrix A°, if N and H are calculated by algo- 
rithm B, then NH = AN where 

^ " ^n-l,(n-l)’V2,(n-2)’ ' ’ ' ^22’^%2’ ‘ ‘ ‘ ^-2,(n-2)’^n-l,(n-l)’ 

and each i' corresponds to the j calculated by algorithm B when r = i - 1. 

Proof — See reference 7. 


Bounds on Rounding Errors 

If the interchanges which need to be made to A are known ahead of time, the 
algorithm for finding N and H from A is the same as algorithm A. 

Up to this point, rounding errors have been ignored. In reality, the computed H 
and N satisfy 


For r = l,2,...,n 

For i = l,2,...,min{r+l,r^ 

n i-1 

= ^ir + I ^ir^kr " I ^ik\r + ®i 

k=2 


hir 


k=r+l 
For i = r+2,r+3,...,n 


ir 


^ir + 


n. 


l,r+l 


11 X 

Z ^ik^^kr - I ^ik^: 

k=r+l k=2 

h. 


kr + ®ir 


r+l,r 


where e^^ is the rounding error made at a given step. 
Let E = then 


NH - AN = E 


( 2 ) 
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Now bounds on E which guarantee that H is very close to a solution of NH = AN 
must be found. Also, 6H must be found such that 


[H ± 6H] c N"^[a ± 


This can be done as follows: Let H' e [H ± 6H] and A' = NH'N then NH' - A'N = 0. 
Subtracting this from equation (2) yields 


N(H - H') = (a - A')n + E 

Therefore 

|A' - a| S |n(H' - H)N'^| + |eN' 1| 
s (n||h - H'l |n~^| + |e||n'1 


Thus 


if 


and 


|a' - a| £ 6A 

|n||6h1|n'^| S — ~ 6A 
R 

IeIIn-^I sisA 


(3) 

(4) 


where R is a large positive number. 

If R is chosen to be very large then j E j must be very small, but 6H gets 
larger as R gets larger. If | E | is very small then a large amount of precision must 
be used in the computations. Since in the algorithms for computing 5H, no attempt is 
made to compute the largest possible 6H satisfying equation (3), there is very little to 
be gained by choosing an extremely large value for R. Usually a satisfactory choice is 
R = 10. Then 6H is taken to be as large as possible and still satisfy equation (3). The 
algorithm used to determine 6H is described in the last section of this chapter. 

Let G = jVjjj where y- > 0 for all i and j. If 

g|n“^|s16A (5) 

then equation (4) can be satisfied by requiring that |e| = G. 
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Let M = N , where M has the same form as N. Therefore, 


G M = 


Ui 

U2 

• • • Un 

^21 

^22 

• • • >^2n 

2^nl 

^n2 

• • • ’^nn 


1 

0 1 
0 I mo 


Since equation (5) must be satisfied, the following equations must be satisfied; 


(i = 1,2, ...,n) 


y* + 

^ir 


3 =r+l 

Equation (6) is satisfied if 


I I < 


(n 5 r > 1; n g i > 1) 


V- . m. — 

^1] jr " R(n - r + 1) 


(i = r,r+l,...,n) 


Therefore, must satisfy 


y.. = min / 

2 ^rSj |R(n - r + l)|mj^ 


From this analysis it is apparent that if is defined to be 


y ■ . = min / 

^ 2grSj|R(n - r + 

= Sail 
R 


(i g 2) 


then, equation (5) is satisfied and |e||m| = ^ for all E such that e^^^ = y^^^. Thus, 
the y^j can serve as bounds on the rounding errors. 
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Effect of Interchanges on Rounding- Error Bounds 

Now, the above assumes that the row and column interchanges that need to be made 
are known beforehand, but these interchanges must be determined by the algorithm; 
and yjj must be calculated before h^j is calculated if i = j + 1 and before n^^ 
is calculated if i > j + 1. 

The following notation is used: b = a @ e means that b is set equal to a value 
which differs from a by no more than e, that is, |b-a| = e; b = a©e means 
that 0 < b - a < e, and b = a © e means that 0<a-b<e; b = a(l © e) is defined 
to be equivalent to b = a © (|a|e) with analogous meanings for b = a(l © e) 
and b = a(l Q e). 

Now, suppose the rounding-error bounds bj^j are determined as in algorithm C. 
Note that it is important that the computed bjj be no larger than the b^ that would be 
calculated if exact arithmetic were used. However, it is not necessary that the b^j be 
calculated extremely accurately since a small change in the rounding- error bound is not 
likely to change significantly the amount of precision required in the corresponding 
computation. Thus, for example, it is required that 


0 0 . 01 ) 

6aii 

that is, b^]^ is required to be less than the exact result of the computation ^ 
but the use of a large amount of precision is not required in the computation. 

Algorithm C 
For i = l,2,...,n 

b.j=£^(i0„.oi) 

For r = 1,2, ...,n 

For i - 1,2, ... ,r 

n 

ir+ I ^ik^kr 
k=r+l 



i-1 

I 

k=2 


i "ik^krl © ^ir 
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If r < n, then 


For i = r+l,r+2,...,n 

n r 

Of = ^ ^ 

k=r+l k=2 

Find j such that 

|»j| =max{|e^^l|, |a,^2|'— |"n|} 
Interchange 

If o'r+i ^ 0, then 

For i = r +l,r+2,...,n 
Interchange 

For i = 1,2 , ...,0 
Interchange 
Interchange (^a- 
Interchange (6a^p6a^+l i) 

For i = 1,2, ....r 
Interchange 

For i =r r+2,r+3,...,n 
^i,r+l ~ ‘^l/°r+l 
^r+l,r ~ '^r+1 
If = 0, then 

For i = r+2,r +3,...,n 
\r+l=0 
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For i = 2,3,...,r 



For i = 1,2, ...,n 



Now it must be determined whether or not the bounds on the rounding errors as 
computed in algorithm C are larger than those given by equation (7). Suppose that 
during execution of algorithm C all calculations for r = r' have been completed. Any 
remaining interchanges affect only rows r'+2,r'+3,...,n of N, A, and 6A, and 
columns r'+2,r'+3,...,n of A and 6A. 

The only elements of M required for the calculation of b-^, (i = 1,2,. ..,n) are 
those in row r'. These elements are not affected by the remaining interchanges. 
Therefore, row r' of M is a row of the inverse of the N obtained after the execu- 
tion of algorithm C is finished. 

For i 5 r' + 1, 6a^^, = 6a. Therefore, for i S r’ + 1, b.^, S y.^,. 

Now, for i > r' + 1, the element in position m ^ after all interchanges are 

completed may be the element which was in position n^, when r = r' during the 

execution of the algorithm. Thus, b., , was used as a bound for the rounding error 

made in the computation of n^ r'+i- This implies that 6a- is equal to the element 

which was in position 6a-, , when r = r'. Therefore, b., , computed during the 

1 , r I , X 

execution of algorithm C satisfies bj^, that is, the bound on the rounding 

error which was made in the calculation of n^ is no larger than the one given by 

equation (7). 

Therefore, if algorithm C is used to compute N and H, 

|e||n'^U-^ (8) 

IX 

where E is the matrix containing the actual rounding errors. 
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An Upper Bound on |m| 

By using variable-precision arithmetic such as SPAR (ref. 8) provides, the ele- 
ments of M can be calculated exactly. However, the bj^j, need not be extremely 
accurate. It is required only that 


°=‘ii 


b.„ i min , ^ 

2^j2r ]R(n - j + l)m 




( 9 ) 


to guarantee that equation (8) is satisfied. Therefore, rather than compute M exactly, 
it is more efficient to calculate an upper bound on |m| by using single-precision 
arithmetic. 

One way of finding an upper bound on |m| is as follows: Let M* be that 

approximation to M which is calculated by 



i-1 


i=p+i 


3|C 


(10) 


where il denotes that the operations are done in single-precision floating-point arith- 
metic. Let Si=nip and = [s^ + ni^j,^p(l + Pk)m*^p^p(l + t7i^)](i + 0j^). Then, 

™ip "" "®i-p some Py., 7?j^, and 0^^ (k = l,2,...,i-p) satisfying |Pj^|, 

< jSgp, and jSgp is a bound on the rounding errors made in single-precision calcu- 
lations (ref. 9, p. 7). Therefore, 

i-p-1 i-1 i -p- 1 

= IT T Tr (1-^ "r) 

i=l 3=P+1 r=j-p 


and 



Therefore 


1-1 

i-1 

i-p-: 

r 

^ z IT 

=1 

j=p+i 

r=j-] 
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Let F = NM* - I where F = If i < j, fij = 0- H i = j, f = 1 - 1 = 0. 

If i = i + 1, f^^ = ^ = 0, since the operation ^ ^ can be per- 

formed exactly. If i > p + 1 

i-1 


'ip=”ip+ 2 "iKp+^i 
j=p+l 


( 12 ) 


Subtracting equation (11) from equation (12) yields 


^ip ’’ip 


i-p-1 

TT (' + »i) 

i=l 


i-1 


I "ii“jp 


j=P+l 


i-p-1 

' - + '=j-p)(‘ + ’’j-p) TT (1 + »r) 

r=i-p 


(13) 


The following lemma is proved on page 65 in reference 10: 
Lemma - K je^| < |3 < 1 and 0 = r S k = n and n/S' = exp / 


- 1 then 


TT (l + ^i) TT + ^i) ^ 1 

i=l i=r+l 


< k/3' 


Applying this lemma to equation (13) yields 

i-1 

|^ip|= Kp|(’ - P - 1)^' +^' ^ 

i=p+i 


n. . 

* 

m. 

1] 

IP 


(i - j + 2) 


where 


Let 


iS' -i| 
n 


n^ 


exp 




1 - 


/3sp/ 


- 1 


i-1 


ip 


nip|(i-p- 1)^’+^' T 
j=P+l 


"ii 


* 

m. 

IP 


(i - j + 2) 


(14) 
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Then 


f. ^ e. 
ip “ iP 


Now 


F = NM* - I = N(M* - M) 
Therefore, for example, 

f — * 

I42 - ™42 " ”^42 


Since |f42| < €43, |m*2 - < e*2- 

Define 643 - ^42- general, for p g q + 2 , 


'pq = V " * "p-q+l ■ ° 

P-1 

-t- V n . (mf - m. ^ 
Zy piwq iqj 

i=q-i-2 


+ m - m 
pq pq 


Therefore 


P -1 


m - m 

pq pq 


= 'pq- I VKq-™lq) 

i=q-i-2 


Therefore 


™pq ■ ”^pq 


^pq 


Vi I ' 
/ n . e- 

Zy pi 1 


i=q-f2 


iq 


Define 


V I I ' 

pq '■pq ■ Z V 

i=qn-2 


I — * 

e„„ = e -1- 


Thus, for p -1- 1 > q. 


m ^ 

< 

* 

m 

pq 


pq 


+ e, 


pq 


( 15 ) 
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Therefore, equation (8) is still satisfied if substitution is made for the part of algo- 
rithm C that governs the calculation of j and as follows: 


For i = 2,3, ...,r 



For i = 1,2, ...,n 



The calculation of m* . . + e' . . requires 3(r - i) multiplications and 
4(r - i) + 1 additions. 

An upper bound on |m| can also be calculated by using expanding -interval arith- 
metic. With this method an interval jm^^ ± 6m^^ is calculated by using the rules of 

expanding -interval arithmetic as described on page 390 in reference 1 or under the 
heading "Rounded Interval Arithmetic” on page 11 in reference 11 such that 
m. e [m.* ± 6m. ” 1 . To find an upper bound on m . . by using this method requires 
4(r - i) multiplications and 4(r - i) + 1 additions. The 6m. calculated by using 
expanding-interval arithmetic satisfies the following inequality: 

^“ip = <p[l + ^sp(i - P - 1)] 
where is defined by equation (15). 
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storage of Arrays 


Although in the above development N, H, and A are thought of as three sepa- 
rate matrices, there is no need to store them in the computer as separate matrices. 
Notice in algorithm C that after a^r is used in the calculation of hj^j, it is never used 
again. Therefore, h^j, can be stored in the same location as a^j,. Also, after a^j. 
is used to calculate n^^ a^^ is not used again. Thus, n^ can be stored in the 

same location as a^^.. If this scheme is used to store N, H, and A, only one set 
of n X n locations is needed to store all three matrices. 


Single-precision arithmetic can be used in the calculation of m*p and Since 

N is calculated by using variable -precision arithmetic, each n^j must be rounded to 
single precision before either m*p or e^q can be calculated. Since the locations 

above the diagonal in the matrix M* = riot used, these locations can be used 

to store the elements of N rounded to single precision. 

Let m?^ = © jSgp) for i > j. With these changes algorithm C becomes 

Algorithm D 


For i = 1,2, ...,n 


^ R 


(l 0 O.Ol) 


For r = 1,2, ...,n 


If r > 1, then 



^ik^k,r-l 


i-1 \ 

I \k-l^kr@^ 

k=2 


^ik^k, 


k=r+l 


,r-l " X 


'kr 


k=2 


If r < n, then 


Find i such that ^ ^|S+l,rl’ K+2,r|’-’ ^nrl} 


17 



If ajjp * 0, then 

For i = 1,2, ...,n 
Interchange 
Interchange 

Interchange (6ajj,6a^^j^j) 
Interchange (6ajj,6a^^^^j) 

^r+l,r ^ ^r+l,r © ^r 


For i = r+2,r+3,...,n 



For i = r+2,r+3,...,n 



For i = 2,3,...,r 
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For i = 1,2, ...,n 


ea„ 


= min > . , , ^ ^ 


(i o o.oi) 


Determination of 6H 

Now that an algorithm for determining N and H has been developed, the 
interval half -width 6H must be determined so that equation (3) is satisfied. Let 


rpqh^.q 

Let 

= [“pq] 

Then 

|m| < M’ 


The problem of finding 6H is solved in two steps. First, a 6Q is found such 

that 

|n|6QsR^6A (16) 

Then 6H is found such that 

6H M’ g 6Q (17) 

Equations (16) and (17) imply that 


|n| 6H |m| 2 |n| 6H M’ ^ |n| 6Q ^ 6A 

R 

Therefore, if equations (16) and (17) are satisfied, then equation (3) is satisfied. 

Let 6Q=j6qjjJ. Because of the form of 6H M', 6q^j = 0 for i s 3. Equa- 
tion (16) is satisfied if 


««ii - ^ 


6q 


^ 6a, 


21 i2 


R 


il 


(i = 2,3, ...,n) 


19 


L 



I I II 


■III 


and 


i 


I l”ik 

k=2 




< R 


R 


6ai- 


(i,j > 1) 


This last equation is satisfied if 


^ik 


R - 1 

5q, • = ^ 

^k] R i _ 1 


(2 S k S i) 


Therefore, if j > 1 and k > 1, define 


SQu = min/^ y- 


"kj R 


igk)|nik|(i - l)j 


Define 


'■Jli = ^ ton 


and 


6q«, = — min/^^ii- 
21 R i£2)K2 


(18) 


(i = l,2,...,n) (19) 


(20) 


With this definition of 6Q, equation (16) is satisfied. 

Notice in equation (18) the quotients ^ k can be calculated beforehand and 

1-1 

stored in the locations reserved for 6a^j. 

If 6H is to satisfy equation (17) then the following equations must be satisfied; 

Shji ^ 6qjj 

ehgj = Sqgj 

and 

n 

I 6hik S 6q. . (j > 1 ) 

k=max^j,i-;^ 
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This last equation is satisfied if 


^hik = 


Therefore, define 


Tk 


aqi^ 

(n - max^i 

- 1} 

■ + ')“ki 

1 X 1 in / 


6q. . ^ 

^1.1 

Kj^k^(n - 

max^ 

j,i - 1} + l)m^.^ 


(max{j,i - ^ k = n) 


(k> 1) 


( 21 ) 


and 


6hii = 6q. j (i = 1 or 2) 


With this definition of 6H, equation (17) is satisfied. 

^^ii 

The quotients — — =5 can be calculated beforehand and stored in the 

n - max^jji -1^ + 1 

location reserved for 6a^j. Once Sq^j is determined, it can be stored in the same 
location as Saj^j if the elements of 6Q are determined in the sequence shown in the 
following algorithm: 


For jr=2,3,...,n 


For k = 2,3, ...,n 


R - 1 ■ ) 

6qi^- = — min< ^ 




The Shj^j^ can be stored in the same location as 6q^j^ if the elements of 6H are 
determined in the sequence shown in the following algorithm: 


For k --= n,n-l,...,2 


For i = 1,2, ...,min{k+l,i^ 


6hik = min 


6q. 


ij 


l<j£k\(n - max{j,i - l) + l)mj^. 


In calculating 6H and 6Q, it is important that the calculated values be no larger 
than the exact values. Putting the above ideas together gives the following algorithm for 
obtaining 6H: 
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Algorithm E 


R' = 


R - 1 
R 


For i - 2,3 ,...,n 

i' = i - 1 
For j=2, 3,...,n 

i 


= 1^)(1 o 0.001) 


6ai, = 


^ O o.ooi) 


li \(n - i + 1) 

6a^j = (R’ 6aj J(l © O.OOl) 


1 6a 


il 


6a«i = R' mia 

' iS2)^|ni2| 

For j=2,3,...,n 


(1 0 0 . 001 ) 


For k = 2,3,...,n 


- 


R' min(^^ 


- max{j,k - 1^ + l) 


(1 0 O.OOl) 


For k = n,n-l, ...,2 


For i = l,2,...,min'{k+l,i^ 


6a,, 


®^ik ^ O o.ooi) 


VKJSkP^kj 


After execution of algorithm E, 6h^j^ is stored in the location initially reserved 
for 6a^^j^. 

This completes the discussion of how to obtain H and 6H. 
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CHAPTER m 


REDUCTION OF AN INTERVAL OF HESSENBERG MATRICES TO AN 
INTERVAL OF COLLEAGUE MATRICES 


Chapter III is devoted to the reduction of an interval of Hessenberg matrices to an 
interval of colleague matrices by a similarity transformation. Where a colleague 
matrix is a matrix of the form 


^1 “1 -^n 

^ ^3 "^n-2 


■ Vl"^2 

1 - ^1 


The interval of colleague matrices is of the form [F ± 6F] where 


0 

0 


6F = • 


0 

0 


6a 


n 



1 


0 ... 0 6aj^ 


A Property of F 

The matrix F has the property that 
n-1 

det(XI - F) = ^ aj^_. Pj(X) + 


( 22 ) 
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where the Pj^(X) are defined by Pj^(X) = det^XI - and 



The proof of this uses the following theorem: 

Theorem 2 - The Pj^(X) satisfy the recurrence relationship 

Pi+l(A) = (X - Pi(X) - 0!i Pi_i(X) (23) 

with Pg = 1 and Pj^ = X - Bj^. 

Proof - Pj.(X) = det(xi - B^) 
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= (x - iS^) det(xi - B^_i) 


+ det 


X - 
-1 


-“1 

^ - h 
-1 


-ac 


-1 X - /Sr-3 -0!r-3 0 

-1 X - ^r-2 0 

-1 -a 


r- 


= (x - Pj.'j det(xi - j det(xi - 


+ det 


X - i3i -ai 

-1 X - |32 
-1 


-«2 


-1 X - |3j._4 0 

-1 X - |3r-3 0 

-1 0 


= (* - <3r) Pr-1<« - “r-1 Pr-2<»> 


since the last term is zero. 

Now it may be proved that equation (22) is valid. 


n-1 

Theo rem 3 - det(XI - F) = ^ a.^_^ P.(X) + Pj^(X) 
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Proof - 


X - -Oil 

-1 ^ ~ ^2 ”“2 
-1 


det(XI - F) = det 


^n -1 

^n -2 


-1 


Vi 

-1 


^2 ■ "n-1 
X - /3n + ai 


-1 ^ - ^2 

-1 >• - ^3 -“3 


= det 


■“n-2 
^ - ^n -1 


X - 


+ det 


-"1 

-1 ^ - ^3 -«3 

-1 A-i 34 


-“n-2 

^ - Pn-1 
-1 


X - jSj -a-^ 

-1 X - ^2 -“2 


+ . . . + (-1)(-Q!jj_j + ag) det 


-1 


^ - V 2 

0 


-O' 


n-2 
-1 _ 


+ (^ - + ^l) 
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Therefore 


det(XI - F) = + a.^- 2^2 + • • • 

+ (”“n-l + ^2)Pn-2 + + ^l)^n-l 

n-1 

^ ^n-i^i " “n-l^n-2 ” ^n)^n-l 

i=0 

n-1 

= > a .p. + P 
n -1 1 n 

i=0 

Because of the recurrence relation in equation (23), the polynomials P^ can be 
made to be the first n polynomials of any desired orthogonal set of polynomials by 
properly selectir^ the ot^ and the jS^. 

Algorithm for Reducing Hessenberg Matrix to Colleague 
Form by Similarity Transformation 

Now in returning to the problem of reducing the interval [H ± 6H] to an interval 
of colleague matrices, it is desirable to determine V and a-j^,a 2 ,...,aj^ such that 


VH = FV 


where V is of the form 



Vll V12 . 


V = 

V22 

V 2 n 



Vnn 


For j < n the ijth element of equation (24) is 
k=i 


(24) 


(25) 
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where 



= V 


n+l,n 


= 0 


and 

V. . =0 (i > j) 

For j = n, equation (24) gives 


I 

1 


Vik\n 


= V 


i-1 


^i^in + “iVl, 


^n-i+l^nn 


(26) 


Define 

istic polynomial of F 


^i,n+l ^n-i+l^nn 

can then be written 


(i = 1,2, ...,n) and 


^n+l,n+l ^nn' 


The character- 


n-1 


det(XI - F) = Pj^(X) + 

Z ^n-i Pi(") 



1=0 


n 



= ^ y 

i4 

Vi^l,n+1 PiW 

(27) 


Thus, the problem of finding is equivalent to finding an additional column 

for the matrix V. 

Equation (25) yields 


^iH-1 = 


3 

k=i 


^kj 


(28) 


3+1,3 


Equation (26) yields 

n 

Vl.n+l = Vj.i „ H- ^ 


(29) 


k=i 


It can be assumed that no h. ^ . is zero since if it is the problem can be treated as 

3+-*-,3 

two smaller problems. 
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Equations (28) and (29) yield the following algorithm for the calculation of V; 


Algorithm F 

Vl,n = 1 


For j = 2,3,...,n+l 



V 1 1 = V 

n+l,n+l nn 


The 


are the rounding errors incurred in the execution of the algorithm. 


Let 



Then 


VH = FV + BK 


(30) 


Obtaining Midpoint of Interval of Colleague Matrices 

Since H is the midpoint of the interval [H ± 6H], it is desirable that F be the 
midpoint of the interval of colleague matrices. How small BK must be and how 
large 6F can be must be determined. 

Let F' be a colleague matrix; then H' = V”^F'V is a Hessenberg matrix, 
where V is the same as in equation (30). Subtracting this from equation (30) gives 

V(H - H') = (F - F’)V + BK 


29 


or 


H - H' = V^(F - F’)V + V"^BK 


V“^Bk|<-^ and 
rt 


IWBK I < which is true if 
R 

(31) 

The proof of this is on page 18 in reference 6. 

A bound on |b- ■ must be known before V.- can be calculated. Therefore, 
from equation (31) it can be seen that the ith column of W must be known before the 
ith row of V can be calculated. There is no direct method of calculating W from V 
so that the ith column of W is available for the calculation of V^. However, W can 
be computed independently of V by the following method: 

If VH = FV then HW = WF. This means that for j £ n - 1 and i = j + 1 

i 

k=max-(l,i-i} 

where w^j = 0 for i > j and w^q = 0. This equation can be solved for w^ 

Thus, W can be determined by the following algorithm: 


It is required that H - H' < 6H, which is true if 


V"^(F - F')v| 6H. 


R 

Define W s v'^. Then it is necessary that 


b. . < — min<; 
D R kSi 


6h 


k] 


(j - k + 2)wj^ih.^i^. 


Wll = l 


For j=2 ,3,...,n 

For i = l,2,...,j 


^ii = 4 , ^ik'^k,j-l - ^i,j-2“i-2 - ^i,j-l^i-l 

k=max-[l,i-ly 


To use the W determined by this algorithm in equation (31) it is required that 
|w| ^ |v”^| but, since the computation of W is done independently of the computation 
of V, this is not necessarily true. If interval arithmetic is used in the calculation of W 
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II 


then it can be guaranteed that |w| exceeds jW'| for any W corresponding exactly 
to any H' e [H ± 6H]; that is, H'W’ = W'F’. 

However, if V is calculated by using algorithm A and if the rounding errors are 
required to satisfy equation (31) with the computed W how can it be guaranteed that 
this V corresponds exactly to some H' e [H ± 6H]? 


The answer to this question is contained in the following theorem which is taken 
from reference 6, page 23: 

Theorem — If W is an upper triangular matrix such that | W| g |W'| for every W 

corresponding to an H’ in the interval [H ± 6H] and if equation (31) is observed in the 

computation of V, that is, if |w||b||k| < then the H' which corresponds to the 

K. 

computed V lies in [H ± 6H]. 

Thus, V”^ corresponds exactly to H' e [H ± 6H]. Therefore, |w|2|v"^|. 

Now, since a bound on |v”^| is known, a bound on jb^^ | can be calculated by using 
equation (31). Thus the following algorithm is obtained for the determination of V: 

Algorithm G 


Vii = l 
wii = 1 

^n+l,n ~ ^ 

For j = 2,3, ...,n+l 


For i = l,2,...,min'(^j,n^ 
If i g n 


[Wy ± aWy] = 


j-1 


4 

k=max{^l,i- 

h-k[v 

'k,i-l 

1} 

- fw, . r, ± 

6W. . r, 

O'* o 

L bi-2 


1 3-2 

■ Ki-i * 






1 

— min< 


6h, 


k,j-l 


R kSi ^(3 - k + l)(|w^| + 6Wi^i)|h.^._j 


(1 G o.oi) 
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'^ii = 


^n+l,n+l ^nn 


/ w \ 

2, V'k.i-l + ’'^1-1, j-l + * “l'^l+l,i-l 

k=l 


\ 




@ ^ 


/ 


that 


Interval arithmetic is used in the calculation of ± Sw-J. 

Obtaining Interval Half -Width 

In the problem of determining how large to make 6F, the constraint on 6F is 


|v"^(F - F’)v| < 6H 
R 


if If - F' < 6F. This is time if 


|W| |F - F’l |V| < 6H 


R 


(32) 


F - FM = 


0 0 
0 0 

0 0 


Therefore 


F - F' V = 


0 0 
0 0 

0 0 


® I ^n " ^n I 

® Kn-1 ■ Vl 


ai-ai 


0 rnnlhn-^ni 

® l^nnl l^n-1 ” ^n-1 


0 Vnn H - ^1 


(33) 


Define ^ ^-i+l^nn ^^i = l^i,n+l ” ^i,n+l 
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Then 


V a . 1 - a 
nn n-i+1 n-i+1 


= V a - 1 -V a'-i 
nn n-i+1 nn n-i+1 


^i,n+l ^i,n+l 


= 6Vi 


( 34 ) 


Up to this point it has been assumed that an interval [F ± 6F] would be obtained, 
but the algorithm determines instead of and it is easier to obtain a 

bound on than on |an_i+l - |. Thus, a bound 6Vj^ is found such that 

if |6Vi| = 6Vi then equation (32) is satisfied. Therefore, n intervals 
[^i,n+l ± (i = l,2,...,n) are found such that if e ± (i = l,2,...,n) 

n 

then the polynomial Pj^(X) + ^ Vl ,n+l ^i- j^(X) is the characteristic poly- 

i=l 

nomial of some matrix H' e [H ± 6H], This is true because, since equations (32) and (31) 

n 

are satisfied, it is implied that Pj^(X) + ^ a^_^ P^(X) is the characteristic polynomial 

i=0 

of some H' e [H ± 6H], but 


'^n+l,„+l P„W 


i=l 


^i-iw = 


nn 


n-1 

i=0 


From equations (33) and (34) it can be seen that 


0 ... 0 

6Vi 


0 ... 0 

6V2 

1 

II 

• 

• 


0 ... 0 

6 Vn 
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Therefore 


|wl If - F'l |v| = 


I 


w 


lil 


w, 


2i 


i=2 


• 0 6V„ 


Thus, if equation (32) is to be satisfied, it must follow that 


n 




hii=V% 


1=] 

This is true if 


5V, |w,d g 


D 1 

< R - 1 jn 


1 I ]i| " R n - j + 1 


that is, if 


6V, S 


< R - 1 




1 R (n - j + l)|w.. 


(j = 1,2, ...,n) 


(i = j,i+l,...,n; j = 1,2,. ..,n) 


0 = l,2,...,i) 


Let 


6V. s min 




^ R Vn - j + 1) |w•^ I 


(35) 


Then = 6Vj^ (i = l,2,...,n) implies that equation (32) is satisfied. 

Qf course, if F and 6F are needed they can very easily be found from 
and 6Vi (i = l,2,...,n). 
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CHAPTER rV 


FACTORING THE CHARACTERISTIC POLYNOMIAL 

In chapter in V^+i^n+l and intervals [Vi^n+l ^ (i = l>2,...,n) were 

obtained such that if Vl e ± 6V^ (i = l,2,...,n) then 

n 

P„(« + I Vl P,.1<X) 
i=l 

is the characteristic polynomial of some matrix H' e [H ± 6H], where the P^ satisfy 
the recurrence relation in equation (23) of chapter HI. 

In chapter IV a method is developed for factoring an interval polynomial into a 
product of interval quadratic polynomials. The notation may be simplified by assuming 
that the interval polynomial to be factored is [P ± 6P] where 

P(x) = a^^PQ + + • • • + ^iPn-1 ^O^n 

and 

aP(x) = Sa^Po + + . . . + 

Also, it is assumed that the are all equal and the jS^ are all equal; that is, 

O'. = ff (i = l,2,...,n-l) 

and 

j3^ = |3 (i = l,2,...,n) 

Given two intervals [C ± 6C] and [D ± 6D], define the product 
[C ± 6C][D ± 6D] = {C'D'IC e [C ± 5C] and D’ e [D ± 6D]}. 

In obtaining a quadratic factor of [P ± 6P], first 

Q = P2 “ 
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and 


= ^O^n-2 + • • • + ^n-2^0 
are found such that 

QT. [p±f] 

Then 6Q and 6T are found such that 

[Q ± 6Q][T ± 6T] C [qt ± 6^ 

that is, 

[Q ± 6Q][T ± 6T] C [P ± 6P] 

The same method applied to [P ± 5 P] can then be applied to [T ± 6T] to obtain 
another quadratic factor. 


Obtaining Interval Midpoints 

Given estimates for s and t, P(x) can be factored in the following manner: 

P(x) = (Pj - sPj - tPo){'=Opn -2 + ■ ■ ■ + + ’’11-2^0) 

+ Vl(Pl-^Po)**>nPo 

To develop a method for factoring P(x) as shown in equation ( 36 ) the following 
lemmas are needed: 

Lemma 1 — If i = 3 ^ 0 then 

i 

P,Pi = ^ aJ->^ P,,2^_j 
k=0 

Proof - For j = 0 and i = 0 , Pg^i ~ ^i’ lemma is true for j = 0 . Now, 

assume 

3 

p.Pl = ^ P,^2^_. 

k=0 
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for i^j if j = For i^m 


PmPi = - »P„.2]Pi 

m-1 

. (X - j Pi^2^.^^1 - a 2 


k=0 

m-1 


m-2 

m-2-k -p 

i+2k-m+l “ M + 2k-m+2 

k=0 


m-1 




-1-k 


k=l 


+ (x - Pj.^^j 


m-1 


i+2k-m 


k=l 


= ^ - /3)Pi^2k-m+l “ “Pi+2k-iJ + ^i-m+l 


k=l 

m-1 




-1-k 


P. 


i+2k-m+2 


+ (x - /3) ^ P 


k=l 

m 


i-m+1 


= I Pl+2k-m ““'\Pl-m+2 + “Pl-m) 

k=2 


Therefore 


m 


P P 
m 




m-k 


P, 


i+2k-m 


k=0 


L emma 2 - P(x) can be factored as shown in equation (36) by using algorithm H. 

Proof - Equation (36) yields 

P(x) = bQP 2 Pj ^_2 + t»iP2Pn_3 ^^2^2^n-4 + • • ■ + V2^0^2 + Vl^l + Vo 
- sbQPjPj^_2 - sbj_PjPj^_3 - ... - sbj^_3P]_Pi - sbj^_2Pi “ ®^n-1^0 


- %PoPn-2 - • • • - %-4 PqP 2 “ t’^n-3PQPl " *^-2^0 
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From lemma 1 


^2^1 = ^i+2 + “^1 + “^^i-2 

and 




Therefore 


P(x) = P„(bo) + P^_i(bi - sbo) + P„_2(bo« + b2 - sbj - tb^) 

+ Pn-sH^l - ^3 - ®^2 - 

+ Pn-4l^0"^ + "(^2 - + ^4 - ®^3 " *^ 2 ] 

+ Pn-sC^l"^ + K^3 ~ ^'" 2 ) + ^5 - ®^4 - 
+ . . . + Pi(bn-5“^ ■*■ “(^n-3 " ®^n-4) *^n-l " ®\-2 ' ^^-3] 

+ PoK^n-4 - "®V3 + " ®^n-l " *^n-2) 

Since the polynomials are linearly independent, their coefficients can be 
equated in the above equation. Thus 

ai = b^ - Sbg 
a-2 =^QOi +b2 - sbj - tbQ 
ag = a^bj - sbQ^ + bg - sbg - tbj 
and for n - 1 g i = 4 

a. = a \_4 + o'(b ._2 - sb._g) + b. - sb._i - tb-.g 


38 



and 


= “\-4 “ \ ‘ ®Vl “ *V2 

Thus the bj^ can be calculated by the following algorithm: 

Algorithm H 

bo = a.Q 

bj = aj + sbg 

bg = a£ + sbj + tbg - abo 

bg = ag + sb2 + tbj - ^(bj - shg) 

For i = 4 ,5,...,n-l 

b. = ai + sb._i + tb._2 - a(b._2 - sb-_3) - a\_^ 

\ = ^K-1 + ' ^\-4 

This completes the proof of lemma 2. 

Let u(s,t) = bj^_ 2 ^ and v(s,t) = b^^. If P 2 - sPj^ - tPQ is to be a quadratic 
factor of P(x) then it is necessary that u(s,t) = 0 and v(s,t) = 0. The problem of 
determining s and t such that this is true is complicated by the fact that it is not 
practical to execute algorithm H exactly. Thus, actually, approximations u* and v* 
to u and v are calculated as in the followii^ algorithm: 

Algorithm I 

bo* = ^0 

— 3,^ + sIDq + 

3|c ^ ^ 

b2 = a2 + sbj + tbg 
bg = ag + sb2 + tbj^ 


- + il^2 

- o(b* - sb*)+ ;//3 
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For i = 4,5,...,n-l 


b* = a. + sb*_i + tb *_2 - «(b *_2 - sb*_ 3 ) - 


= ^K-1 + ^^n-2 + ^K-3 - ^\-4 + ’^n 

u*(s,t) =b*_j 
v*(s,t) =b* 

where the i{/^ are rounding errors. Let a! = a^ + Then 


^ 0 Pn + ^lPn-l + - • - + ^ 


;Pq = (P2 - "Pi - tPo)N*Pn-2 + 
+^*-i(pi-"Po)-^^;po 


^VsPi-^b 


This can be rewritten as 


P’(x) = (?2 - sPj - tPo)(bQP„_2 + . . . + b*_3Pi + b^_2Po 


where 


P’(x) = agP^ + a'jP^_^ + + a;_ 2 P 2 + Pi(^n-1 ' '"n-l) + PqK " K + "C 


Then 


if 


P' e 


P ± 


a. -ai 


R 


6 a^ 


(i = l,2,...,n-2) 


6 a 


n-1 


Vl'^n-l"Vl|= R 


and 


!a' - b* + sb* - a I ^ 
i^n °n ^ "°n-l ^n| " r 
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3 * 



These equations are satisfied if 


6ai 




IVil ■ 


6a, 


n-1 


2R 


, 5a„ 

I 


b 


n-1 1 


< ^^n-1 


2R 


(i= l,2,...,n-2) 


and 




6a„ 
< n 

" ~2R 


The above equations can be checked to determine if improved estimates s and t 
are needed. New estimates Sj^ and tj,^ can be calculated by using Newton' s method; 
that is, 


where 


'®N 


's' 


'As 


= 


+ 




t 


At 


as 

au 

at 


As 


U 

3v 

_3s 

3v 

9t_ 

s.t 

'> 

L . 


_t_ 


Approximations and to — and 

as as as 

following algorithm: 


— , respectively, can be computed by the 

OS 


Algorithm J 



dj = bj + sdQ + rjj 

d* = h*+ Sd* + tdj - Q! (d* - b*) + 772 
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For i = 3,4,...,n-2 


df = bi + + tdf_2 + «(bf_2 - df_2 + sdf_3) - « d*,^ + r?. 


^n-1 = '"n-l + K-2 + K-3 + "^n-3 + "^^4 " ^\-5 


9 u 


^ = V2 


— =d: 1 

8s 


The ri^ are rounding errors. If each = 0 and if h* = then 


and 


dv* _ ^ 
8s 8s’ 


8u* ^ ^ 
8s 8s 


Algorithm K can be used to compute approximations to and — as follows: 

8 t ot 


Algorithm K 


c!a = 0 

C*1 = 0 
Co* = b* 

C*=b* + sC*+q 
For i - 2,3, ■■.,n-3 


2^* 


C* = b* + SC*_1 + tc*_2 - »(C*_2 - sc'.j) - + 5 j 


= b!_, + sc*_3 + tc:_^ + asc‘_^ - a^C* 


■'n-2 ■ ''n-2 


n-4 


"n-5 


"n-6 


8 u _ p* 

8t ^ n-3 


8v _ p* 

8t '"n-2 


* 

The are the rounding errors and and 


8v 


8u 


8v 


are approximations to — and — , 
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Then the As and At satisfy 


8u* 

0s 

0V* 

0S 


0u* 

0t 

0V* 

at 



As - Pj 

1 

11 

u* 

s,t 

At-Pa 


v*_ 

l_ — 




s,t 


where Pj and Pg take into account the rounding errors made in determining As 
and At. Thus, the new estimates for s and t are Sj^^ = s + As and tjq^ = t + At. 

Successive estimates for s and t can be calculated in this manner, but the 
process does not cause convergence unless the rounding errors are controlled. The 
following method is chosen for controlling rounding errors: As the (i + l)th 


estimate 


®i+l 

^i+1 


is being calculated, each of the rounding errors [Vj. (r = l,2,...,n). 


7]^ (r = l,2,...,n-l), (r = 1,2, ...,n-2), and Pi,P^ is required to be in absolute 

2 ^ 

, where X = . — and 


value less than X 


^i ■ ®i-l 
ti - ti-1 


So 

L^oj 


^0 


is the initial estimate. 


Thus the algorithm for obtaining the interval midpoints is: 
Algorithm L 


X = 


max{|so|, |tol} 

For j = 0,1,..., large 


(l 0 O.Ol) 


S = S; 


* 


^0=^0 


b* = (aj + sb*) 0 I 

bg = (cL2 + sbj + tbg - Q-bg^ 0 ^ 

ag + sbg + tbj - a(^bj - sbg j 


b* - 
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For i = 4,5,...,n-l 


b* = [aj + sb*_i + tb*_2 - ^(bf.g - sb*_g) - q-\*_4] © 4 

= k + K-1 ^ K-2 + "<-3 - “^^ 4 ) @ ^ 


“ K-1 


6a. 


n-1 

2R 


6a„ 




|6a-| 6a<] o 1 

. J 1 2 n-2 n-1 nv,., 

R R 2R 2Rj 

convergence has been obtained. Proceed to the 
determination of the interval half-widths. 


d!i = 0 

^0 = ^0 


dj = (b* + sd*) © I 
For i = 3,4,...,n-2 


© ^ 


d* = [b* + sd*_j + td*_2 + «(b*_2 - d*_2 + sd*_3) - 
<-l = K-1 + ®<-2 + K-3 + ^K-3 + "<-4 - “^<- 5 ) © ^ 


C_*2 = 0 
C_i = 0 


r* - h* 

Uq - Dq 


cj = (bj* + scj) @ J 
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For i = 2,3,...,n-3 


C* = [b* sC*_i + tC*_2 - «(C*_2 - SC*_3) - a^cl^ © ^ 

^n-2 = fe-2 + ®^n-3 + + "®^n-5 " "^^n-e) © ^ 


bn-lCn-2 “ bnC3 
^ ‘^11-2^2 ■ %-l^n-3y 


As = 


At = 


K <-2 - bn-lCl 

<-2<-2 - <-K-2J 


= x(min{|As|, |At|}) 


Vi " "r 


Vl = V 


Let 





X-^ 


s 

x2_ 


_t_ 


L=X X, 


X 


i-1 


2 


and 


Let 


f(x) = 



f*(x) 


u*(s,t) 

y*(s,t)_ 


8u &u 

as at 

ay 8v 

as at 
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and 


J* = 


9u* 

9u* 

9s 

9t 

9v* 

9v* 

9s 

9t 


Then, theorem 4 may be stated as follows; 

Theorem 4 - The convergence criteria of algorithm L are satisfied in a finite number 
of iterations if the following assumptions are satisfied: 

(1) J~Vo) 

' g b 


( 2 ) 

(3) X 


Xi - Xq 


k=l 


d\(x) 


9x^3x^ 


< £ 
" 2 


for all X in ||x - Xq|| S 2b where x^ = s, x^ = t, f ^ = u, and f2 = v 


(4) h < i where h = abc 

2 

1 1 '^ 

(5) XabK < i 


where 


and 


2 5 + 8h 


K = sup k(x) 
||x-Xo||S2b 


k(x) = i + 6f(x) + 6j(x) 

a 1 - 2h 


where 6f and 6J are positive continuous functions of x with the property that 
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and, for each i, 


Proof 


bQ,bi, 


Let 


and 


Let 


Then 


f *(>4) - '(>'1)11 - « “(>4) 


(6) Xab^ 


sup 6J (x) 

jx-xo!1^2b 


<i 

8 


(7) ^0 ^ 4Xb^ 

— See appendix A. 

Obtaining Interval Half -Widths 

From now on the bQ,b^,...,b*_2 calculated by algorithm L are denoted by 
•••’^n-2- 

The interval midpoints s and t (bQ,bj,...,bj^_ 2 ) satisfy 

^O^n ^ ^'l^n-1 + • • • + ^n-2^2 ^l(^n-l " ’^n-l) ^o(^n " '’n ^’^n-l) 

= (P2 - sPi - tPo)(boP^_2 + . . . + b^_gPi + b^_2Po) 


a. =a. 


(i = l,2,...,n-2) 


a" 1 = a’ - b , 
n- 1 n- 1 n- 1 


^n = ^n - ^n + ®Vl 


n 


P" = »0Pn ■> 2 "I'P! 


i=l 


P e 


P ± 


R 
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Consider 


Now, if 


P(x) = |P2 - (s + 6s)Pi - (t + 6t^[bj^_2 + 6V2 +(^n -3 + ^Vs)^! 
+ •••■*• (bi + ^^l)Pn -3 + ^O^n-2] 

P(x) = + 6a|^ + + 6a|^_i)Pi + • • • + (a'j + 5a'j)Pj^_j + apP^^ 


then expanding the first equation for P(x) and equating coefficients of the Pi(x) as 
was done in lemma 2 yield 

ao = bQ 

a'j + 6 a j = bj + 6bj^-(s + 6 s) bQ 

a 2 + 6 a 2 = bgO! + b 2 + 6 b 2 - (s + 6 s)(bj + 6 b - (t + 6 t)bQ 

ag + 6ag = Q!jbj + 6b j - (s + 6s)boj + bg + 6bg 

- (s + 6 s)^b 2 + 6 b 2 ^ - (t + 6 t)^bj^ + 6 b^^^ 

and for n - 2 S i g 4 

a!' + 6aV = a^(b^_4 + Sb-^^) + Q!|bj_2 + Sb- g - (s + 6s)(b._g + 6b-_3) 

+ (b- + 6 b.) - (s + 6 s)(bj _4 + 6b. _j) - (t + 6 t)(b ^_2 + Sbj^_ 2 ) 

^n -1 + K -1 = (Vs + V-s)“^ + "l^n-3 + ^V3 ‘ + ^®KV4 + ^V4) 

- (s + Ss)(bj ^_2 + 6 bj^_ 2 ) - (t + 6 t)(bj ^_3 + 6 bj^_ 3 ) 

^n + K = A\-4 + ^- 4 ) - + ^«)(V3 + ^^n- 3 ) - + ^^H'^n-2 + V- 2 ) 

Using the fact that 

P"(x) = (P 2 - sP, - tPo)(b„P „_2 + . . . + b^.jPi + b„_ 3 P„) 
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yields 


and for n - 


Now i 
such that if 

P(x) e |p" d 


For n - 2 


6aj = 6bj - 6s bg 




eag = 6 b 2 - s 6 bj - 6 s(bj + 6 bj^ - 6 t bQ 
6 ag = Q!( 6 bj - 6 s bg) + 6 bg - s 6 b 2 - 6 s(b 2 + 6 b 2 ^ 

- t 6 bj - 6 t(bj^ + 6 bj) 

2 5 i > 4 

6aj’ = Q!^6bj_^ + Q!j6b^_2 ~ ® ®^i ” ® ®^i-l 

- 6 s(b._j + 6 b. _j) - t 6 b . _2 - 6 t(b ._2 + 6 b. _ 2 ) 


( 37 ) 




S^n-3 - ® ®^n-4 ” ^K^n-4 + ®^n-4) 


- ® ^^n-2 - ^K^n-2 + ®V2) ’ ‘ ^^n-3 ' ®K^n-3 + ^K-i) 

6 a” = 6 bj ^_4 - O' 6s(b^_g + 6b^_g) - sa 6b^^_g - t 6 b ^_2 - 6t(bj^_g + 6 b^^_ 2 ) 

t is desirable to find bounds 6s on 6s, 6t on 6t, and 6b^ on 6b^ 

|6s| = 6s, |6t|= 6t, and |6bi| = | 6b-[ | (i = l,2,...,n-2) then 

R - 1 


R 


6P 


. This is true if 


Jsrt ^ ^ <. “ 1 C 


(i = l,2,...,n) 


(38) 


i = 4, where M > 1, equations (38) are satisfied if 


6b, 


i-4 


R^6a. 

5M R 1 


s 6b, 


M - 1 R - 1 


6a, 


i-3l “ 5M ■ R “"1 

M - 1 R - 1 


(|q!| + |t|) 6b._2 


5M R 
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ci-l-eM -1 R- 1 r„ 


l«IN I Vs + ®Vs| = 


-L 5 ^ 6 a. 

3M R 1 


|6s| bj_i + 6b._j 


< 1 R - 1 6a 
" 3M R ® 1 


! 6 tl|bj _2 + 6 bj. 2 | 


Treating the equations for 6aj, Sag, Sag, Sa^^j, and Sa^^ in: 
shows that equations ( 38 ) are valid for |Sb^| s sb-, |Ss| S Ss, and |st| 
Sbj^,...,^^ 2> Ss, and St defined in the following manner: 

Sbo = 0 (i £ 1) 

^ 5M R I ^ |s| |t| + |o!| |s| q|2 


where Saj = 0 for j > n. 


St = ^ min 


Sa. 


3M R n^i^2 1 |bj_g + Sb._2| 


Ss s -1- — 1 mW min 


Sa^ 


3M R 


min 


Saj 


n^i^nVj +SVij n5i53|^|o!|( b._g| 4 


Let 


Ql - Pg - SPj - tPg 
SQj = Ss Pj + St Pq 

= ^oV -2 + ^lV -3 + • • • + V 2 V 


1 similar manner 
= St with 


(39) 


(40) 
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and 


Then 


6Ri = 6bl Pn_3 + P„_4 + . . . + 6b^_2 Pq 


|Qi ± 6QJ[Rj ± 6Rjj c [p ± 6p] 


The same methods can be applied to jkj ± 6 rJ to calculate [Qg ± SQg] 
and ± 6R2^ such that 

[Q2 ± 6Q2][R2 ± 6R2 ] C [Rj ± 6 rJ 

This process can be continued until 


± 6 qJ 1^2 ± • ■ • 

Q 

n' 

± 6Q 

fil 

R 

n’ 

± 6R 

n' 

C [p ± 6p] 



2 


_2j 


2 


2 _ 



where Ri-j^-, is either a constant or a polynomial of degree one depending on whether n 
2 

is even or odd. 

The roots of any Qj e [Q^ ± are eigenvalues of some matrix A' e [a ± 6 a]. 
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CHAPTER V 


CONCLUDING REMARKS 


In chapter II algorithms are given for reducing an interval [A ± 6 A] of real 
matrices to an interval [H ± 6H] of Hessenberg matrices such that each element 
of [H ± 6H] is similar to some matrix belonging to the interval [A ± 6A] . The algo- 
rithms of chapter II can be used in developing contracting- interval programs for other 
processes which require reduction of a matrix to Hessenberg form. 


In chapter III algorithms are given for reducing an interval of Hessenberg matrices 
to an interval of colleague matrices by a similarity transformation. As is shown in 
chapter III this interval of colleague matrices is equivalent to an interval of character- 
istic polynomials [P ± 6P], each element of which is the characteristic polynomial of 
some H' e [H ± 6H]. 


In chapter IV algorithms are given for obtaining interval quadratic poly- 


nomials ± 6Q]j, |Q2 • 5Q2], • • •, 
|qj ± sQjj |q£ ± 5Q2J . . . 




± 6Q 


such that 


Q 

n’ 

± 6Q 

'nl 

R 

n' 

± 6R 

n"] 


2. 


2j 


2. 


2J 


Q [P ± 6P] 


where 

odd. 



is a constant if n is even or R is a polynomial of degree one if n is 


Now if the Hessenberg matrix H that is obtained from chapter II has any zeros 
along the subdiagonal then it is partitioned into Hessenberg matrices Hj, H2, . . ., H^ 


with 



order (H^^) and such that no Hj has a zero along the subdiagonal. 


Each Hj^ must be reduced to colleague form by the methods of chapter HI, and then 
the characteristic polynomial factored by the methods of chapter IV. 


Combining the methods of chapters II to IV yields quadratic interval poly- 


nomials 


Qy ± aQ,jj 


i = 1 , 2 ,..., 


; i = l,2,...,zj and interval polynomials j^R^^ ± SrJ 

(i = where Rj is a polynomial of degree one if order (Hj) is odd and Rj is 

a constant if order (H^) is even. Then if 





i = l,2,. 
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and 


then 


e [Rj ± 6Ri] 


/ 


I 

i=l 


R. 


•• TT Q'-- 

\ i=i / 


(i = 1,2,...,Z) 


(42) 


is the characteristic polynomial of some matrix A' e [A ± 6A]. 

Thus, any root of equation (42) is an eigenvalue of some matrix A' e [A ± 6A]. 
Therefore, since the roots of any Q^j are trivial to determine, a method is developed 
for obtaining eigenvalues which are exact for a matrix differing by less than a specified 
amount from the matrix A. 

Appendix B contains an example of this method applied to a simple 3x3 matrix. 

It is essential that variable-precision arithmetic be used when programing the 
algorithms given in this report. It cannot be assured that the program is a contracting- 
interval program unless the required amount of precision is used, and this required 
amount of precision may be more than is obtained with either single- or double-precision 
arithmetic. 

Variable-precision arithmetic is presently not available at very many computing 
installations since the hardware on most present-generation computers was not designed 
in a way which would facilitate implementation of variable -precision arithmetic. How- 
ever, software implementations of variable-precision arithmetic such as SPAR (ref. 8) 
have been developed. As more efficient implementations of variable-precision arith- 
metic are developed, the advantages of contracting- interval programs will be more dis- 
cernible. Even with efficient implementations of variable-precision arithmetic it will 
probably talce more computer time to execute an algorithm by using variable-precision 
arithmetic than by using single-precision arithmetic. However, this should not lead to 
the dismissal of the concept of contracting-interval programs. The accuracy of computed 
results must be determined by some means if these results are to be of any use, and it 
may be more efficient to use a few more minutes of computer time to calculate results of 
known accuracy than to estimate the accuracy by some other means. Conventional 
methods of determining accuracy may require a comparison of the results of several 
computer runs or a study of the problem by a numerical analyst. 
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It has been the accepted pattern to have computers assume more and more duties 
previously done by people, so the concept of computing numbers of known accuracy 
should be the natural thing to do when the computers are designed with this in mind and 
when the methods become available. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., July 29, 1971. 
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APPENDIX A 


PROOF OF THEOREM 4 

The proof of theorem 4 makes use of a theorem which gives sufficient conditions 
for the convergence of Newton's method in R^, where the function and its Jacobian are 
known only approximately. This theorem is a modification of a theorem given on 
page 115 in reference 12. 

Suppose it is desired to find a root of f(x) = 0, where 



fl(x) 


xl' 

f(x) = 

f2^ 

and X = 

x2 


fn(x) 


xn 


By using Newton's method successive approximations would be generated by 
solving J(xJ (^i+1 " ^i) ~ "^(^i) ^i+1’ '^(^i) Jacobian evaluated 

at x^. In most cases neither J nor f can be evaluated exactly at the point x^. 

Let J. (x.) and f^ (xj) be the approximations to J(x^) and f(xj) that are calculated. 

Therefore, actually, 

J|i(^i) (^i+1 - ^i) = -f?i(^i) 


is solved for x^^^j except that this system of equations cannot be solved exactly for 
^i+1' Therefore, x^^^^ satisfies 

J?i(^i)(Vi-^i-^i) = -f?i(^i) 


where is the error made in the calculation. 

Now, if the approximations J^ (x^^) and f|.(xj^) 

or if Pj does not approach 0 as i approaches <» 
not be expected to converge to a root of f(x) = 0. 


do not approach J(x^) and f(x^) 
then the sequence of iterates can- 
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APPENDIX A — Continued 


However, suppose it can be guaranteed that 


and 


J|i(^i)- J(^i)|L= 

lhlL=«i 


(A2) 

(A3) 

(A4) 


where 6J(x) and 6f(x) are positive continuous functions of x, 




li = X Xi-x._i 


is a prescribed value, and X is a constant having the same dimensions as 1/x. 

The theorem gives sufficient conditions for the convergence of the iterates deter- 
mined by equation (Al) where equations (A2) to (A4) are satisfied at each step of the 
procedure. 

Theorem 5 - Let successive approximations to a root of f(x) =0 be generated by 
solving equation (Al) for x^^^j at each step, where equations (A2) to (A4) are satisfied 
and Iq g 4Xb^. 

All norms used are °°-norms. The following assumptions are made; 

(1) If Xq is the initial iterate, then jjj”^(xQ)|j = a 

(2) ||xi-XojUb 


(3) Let the components of f(x) have continuous second derivatives which satisfy 


XX 

I 

k=l 


a2fi(x) 


ax^8x^ 


= r- for all X in x - Xnfl = 2b 


( 4 ) 

( 5 ) 


Let h < i where h = abc 
2 


Let f^(x) be bounded for x - Xq S 2b and ^ g ^q, and let 


F s sup |jffc(x)|| 
||x-xo||S2b" ^ 

mo 
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APPENDIX A — Continued 


Let 


and 


k(x) s i + 6f(x) + 6J(x) 

K = sup k(x) 

||x-Xo||g2b 


Then it is required that 


XabK < - 
2 



5 + 8h 


( 6 ) 


Xab^ 


sup 6 J (x) 

lx-xo||S2b 



If all these assumptions are satisfied then the iterates are uniquely defined and 
II X - xqII = 2b for each v; and the iterates converge to some vector, say lim x = a, 

V-^oo 

for which 


f(o!) = 0 

and 

IK-“II = P 

Proof - See reference 7, page 11. 


i-1 

Lemma 3 — Let = y^ + ^ ajj(s,t) 

^ i=i 

i-1 

w! =y:+ ^ Qfi.(s,t) w! 

i=i 

and Zj = - wj^ (i = l,2,...,n) where each a^j(s,t) is a continuous function of s 

and t. If Cj^(s,t) is a positive continuous function of s and t for each i and 
if |yj - y|| < ^ Cj(s,t) (i = l,2,...,n) then jz;^| < § Bj{s,t) for some positive continuous 
function B^(s,t). 
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APPENDIX A — Continued 


Proof — For i = 1 define = Cj^ and the result is obvious. 
Assume |Zj^| < (i = 2,3,...,k). Then 


‘k+1 


yfcH.1 - yLi * I “k+i,f j 

i=i 

k 


= K+1 - yk+il 


i=i 


k 

Z |"k+l,j|h| <^<^k+l+ ^ |"k+l,j|^®j 

j=l 


Define Bj^^j(s,t) b Cj^+i(s,t) + ^ .(s,t) B.(s,t). Then 

i=i 


and is a continuous function of s and t. Let x = 


. Then 


and 


f(x) = 


u(s,t) 

y(s,t)_ 


Let 


f*(x) 


u*(s,t) 

y*(s,t)_ 


9u 

9u 

9s 

9t 

9v 

9v 

9s 

9t 


and 


9u* 

9u* 

9s 

9t 

9v* 

9v* 

9s 

9t 


Theorem 6 — Let f*(x) be calculated by algorithm I where | (i = l,2,...,n). 

Then there exists a positive continuous function 6f(x) such that ||f*(x) - f(x)|| < ^ 6f(x) 
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APPENDIX A - Continued 


for all X. Let J*(x) be calculated by algorithms J and K where |77^j < | 

(i = l,2,...,n-l) and |zj| < | (i = l,2,...,n-2). Then there exists a positive continuous 

function 6J(x) such that ||j*(x) - J(x)jj < | 6j(x) for all x. 

Proof - Since the multipliers of the bj and b? in algorithms H and I are continuous 
functions of s and t and, since 


ai - (aj + ^ 


(i = l,2,.,.,n) 


lemma 3 can be applied to obtain |b^ - b*| < ^ Bj(s,t) for some positive continuous 
fimction B^(s,t). 


Let 


6f(x) = 


B„_l(s,t) 

Bn(s,t) 


then 6f(x) is a continuous function of x and 


|f*(x) - f(x)|| = 


u*(s,t) - u(s,t) 
_v*(s,t) - v(s,t)_ 




< 4 Sf(x) 


Now 




bf - b?| + |T?i| S ^[Bi(s,t) 



Let d^ be the result of algorithm J when each = 0. Then by lemma 3 there exist 
positive continuous functions D^(s,t) such that 


dj - df < ? Di(s,t) 


Therefore 


and 


8u 9u* 
9s 9s 


8v 9v* 
8s 9s 


< I Dj,_2(s,t) 

< ?Dj^_i(s,t) 
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II II III 


APPENDIX A — Continued 

Likewise, if Cj is the result of algorithm K when each = 0 

|ci - cj'l < I Ci(s,t) (i = l,2,...,n-2) 

for some positive continuous functions Cj(s,t) (i = l,2,...,n-2). Therefore 

< I Cj^_3(s,t) 

< I Cj^_2(s,t) 


8u 9u* 

at ‘ at 

and 

dv dv* 

W ~ ^ 


Define 


6J(x) 


Dj^_2(s,t) 

C„_3(s,t)1 


C„_2(s,t) 


Then, 6J(x) is a positive continuous function of x and 
|j*(x) - J(x)| < ^ 6J(x) 


Lemma 4 - In any bounded region ||x - Xq[| < r, if 4 < ^0 where is some positive 

value, then there exists a constant F such that j|f*(x)|| < F for each x in this region. 

Proof - |)f*(x)|l < l|f(x)|j + 6f(x). The result follows from the fact that both f(x) 

and 6f(x) are continuous functions. 

Theorem 7 - The convergence criteria of algorithm L are satisfied in a finite number of 
iterations if the following assumptions are satisfied: 


( 1 ) 

( 2 ) 

(3) 


J-l(xo)|i §a 


|xi - xo| S b 


I 

k=l 


a2f.(x) 

ax^ax^ 


< 


c 

2 


for all X in 


1 ^ 

X - Xq[ ^ 2b where x = s, x = t. 


f j = u, and f 2 = V 
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( 4 ) 


h < — where h = abc 
2 


(5) 


1 5-'' 

XabK < i 


2 5 + 8 h 


where 


and 


K = sup k(x) 

l|x-Xoll^2b 


k(x) = i + 6 f(x) + 6 J(x) 

(6) Xab^ / sup 6J (x)\ < i 

y|x-Xo||S 2 b j 8 

( 7 ) £ 4 Xb^ 


Proof — In iteration k of algorithm L the rounding errors made in calculating the b*, 
Cp and d* are all less than Thus, theorem 6 can be applied to get 


and 


|f*(^k) - f(^k) 


< ^k ®f(^k) 


|j*(Xk) - J(Xk)||< ^J(^k) 


The rounding errors made in calculating As and At are less than Thus, 

all the requirements of theorem 5 are met. Therefore, the successive iterates Xj con- 
verge to some a for which f(o!) =0. Also, since x.; — a;, ^ 0. 

J J 

f 6 ai 6 ao 6 an _9 6 a„ i 6 a„"^ 

Let X = min|^, There exists Nj such 

that j > Nj implies that < X. 

Now 0 = f(a) = f/lim xA = lim f('x.V Therefore, there exists No such that 
\j^°° 7 ^ 

for j > N 2 
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I 

''h’‘i)| ^ 4R 


U(Sj,tj)| 


< min 


^^n-1 ^^n 

. 2R ’ |s|4R 


These two equations imply that for j > N£ 




/ 

v(Sj,tj)-Su(Sj,tj) 


Therefore, for j > max^N^jNg^ the convergence criteria of algorithm L are met. 
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APPENDIX B 


EXAMPLE 


In appendix B, the algorithms of this report are applied to a 3 x 3 matrix. It is 
hoped that this will aid the reader in understanding which algorithms to use, the sequence 
in which they are applied, and the result of each algorithm. Let 



5 

-6 

-6 


0.1 

0.1 

0.1 

A = 

-1 

4 

2 

6A = 

0.1 

0.1 

0.1 


3 

-6 

-4 


0.1 

0.1 

0.1 


R = 10 and each 0!p/3j = 0. Use of algorithm D yields 


5 

-4.02 

-6 

0.3 

-2.02 

-6 

0.33 

0.01 

2.02 


and use of algorithm E yields 


6A = 


0.09 

0.045 

0.09 


0.09 

0.045 

0.09 


_ 

0.0225 

0.045 


Now 

let the 

upper Hessenberg portion of A be renamed 


H = 


-4.02 

-6 

-2.02 

-6 

0.01 

2.02 


and let 6H = 6A. Then, algorithm G yields 


1 -1.667 65.266 


V = 


0.3333 -99.3734 


33.33 
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The fourth column of V is 

-135.8393 
268.0000 
-166.7000 
33.33 

Equation (35) yields 

= 0.027 
eVg = 0.0054 
6V3 = 0.0015 

Thus, the interval polynomial is 

33.33X^ + [-166.7 0 0.0015]x^ + [268.0 0 0.0054]x + [-135.8393 0 0.027] 

Now algorithm L is used to factor this polynomial. If sq = 3, tQ = -2, and = lO”"^, 
then after the first iteration Sj = 3.0022, tj = -2.0341, and = 10“^. After two 
iterations S2 = 3.00211, t2 = -2.03841, and ^2 ~ 10"®. During the third iteration there 
are computed 

bg = 33.33 
b] = -66.63967370 
bj = 0.00016389 
bg = 0.00016928 

The convergence criteria are met. 

Now, by using equations (39) to (41) and with M = 3, 6bj, 6t, and 6s are com- 
puted as 

6bi = 0.00018 

^ = 0.00004051 

and 

^ = 0.0000081 



64 



APPENDIX B - Continued 


Thus 


(x^ - [3.00211 0 0.000008l]x - [-2.03841 0 0.00004051]) (33. 33X + [66.6396737 0 0.0001^) 
c 33.33X^ + [166.7 0 O.OOlsjx^ + [268.0 0 0.0054]x + [135.8393 0 0.027] 
Note that, since 


H = 


-4.02 

-6 

-2.02 

-6 

0.01 

2.02 


and 


6 H = 


0.09 

0.09 


0.045 

0.09 

0.045 

0.09 

0.0225 

0.045 


if the 0,01 is replaced with 0 in H and the 0.0225 replaced with 0.0125 in 6 H, 
since [O 0 0.012^ c [O.Ol 0 0.0225], then 


H = 


5 -4.02 -6 

3 -2.02 -6 

2.02 


and 



0.09 

0.045 

0.09 

6 H = 

0.09 

0.045 

0.09 



0.0125 

0.045 


Then H can be partitioned into submatrices which can be treated separately. Let 



-4.02 


- 2.02 


6Hj = 


0.09 

0.09 


0.045 

0.045 


H 2 = 2.02 and 6 H 2 = 0.045. Thus, for \ e [2.02 ± 0.045], X is an eigenvalue of some 
matrix in [A ± 6 A]. 
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Using algorithm G on H]_ yields 



1.667 

0.3333 


The third column of V is 

0.652 

-0.9937 

0.3333 


6Vi = 0.02025 

and 

^2 = 0-0081 


Therefore, the characteristic polynomial of [h^ ± 6HjJ is 

0.3333X^ + fo.9937 ± 0.008l]x + [0.652 ± O.O2025] 


Thus 


(x - [2.02 ± 0.04^)(0.3333X^ + [-0.9937 ± 0.008l]x + [0.652 ± 0.02025]) 


is an interval polynomial, each element of which is the characteristic polynomial of 
some A' e [A ± 6A]. 
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